/**
 * visHierarchicalClustering.h
 * 7/4/2009.
 * linux(ubuntu8.10), g++4.3.2
 * Author : luxinlai
 */

#include "visHierarchicalClustering.h"

//#define USE_MIN_DIS_DIRECTLY

//#define USE_MIN_DIS

//#define USE_AGGRESSIVE
//#define USE_EU_DIS
//#define POST_FILTER

//#define USE_EU_NO_COMBINE

//#define Force_Init_Seeds_For_Topo_Clutering


bool *longLine;
namespace visualization
{
void visHierarchicalClustering::initialize()
{
	m_ObjNumber = 0;
	const double smallSize = 5;
	longLine = new bool[getStreamlinesNumber()];
	for(size_t i = 0; i < getStreamlinesNumber(); i++)
	{
		if(smallSize > getStreamlineSize(i))
		{
			longLine[i] = false;
			continue;
		}
		longLine[i] = true;
		m_ObjNumber++;
	}
	cout << "Actual streamline number	:	" << m_ObjNumber << endl;
	m_SimilarityMatrix = new double*[m_ObjNumber];
	for(size_t i = 0; i < m_ObjNumber; i++)	
		m_SimilarityMatrix[i] = new double[m_ObjNumber];

	computeLength();

	m_Centers = new vector<Point>[m_ObjNumber];
	size_t indx = 0;
	cout << "Initialize centers..." << endl;
	for(size_t line = 0; line < getStreamlinesNumber(); line++)
	{
		if(!longLine[line])continue; // not long line
		for(size_t sz = 0; sz <= m_seg; sz++)
		//for(size_t sz = 0; sz < getStreamlineSize(line); sz++)
		{
			Point pnt;
			getPointOnStreamline(line, m_seg, sz, pnt);
			//pnt = getPointOnStreamline(line, sz);
			m_Centers[indx].push_back( pnt );
		}
		indx++;
	}
	m_PartsNow = m_ObjNumber;
	m_IsCenter = new bool[m_ObjNumber];
	for(size_t i = 0; i < m_ObjNumber; i++){ m_IsCenter[i] = true;}

#ifdef USE_MIN_DIS
	buildMinDisTree();
	cout << "Computing minimum distance..." << endl;
	for(size_t i = 0; i < m_ObjNumber; i++)
	{
		for(size_t j = 0; j < m_ObjNumber; j++)
		{
			m_SimilarityMatrix[i][j] = minDisSim(i, j);
		}
	}
#endif


//	m_Length = new double[m_ObjNumber];
//	size_t index = 0;
//	for(size_t i = 0; i < getStreamlinesNumber(); i++)
//	{	if(longLine[i])m_Length[index++] = m_streamlineLength[i][ m_streamlineLength[i].size() - 1 ];}

}

void visHierarchicalClustering::partition(size_t &part, double predist, double distance)
{
//Eu_Filter(part, predist);
//return;
filter(part, predist);
//return;

	size_t indx1 = 0;
	cout << "Compute initial similarity..." << endl;
#ifdef Force_Init_Seeds_For_Topo_Clutering
	for(size_t r = 1; r < m_ObjNumber; r++)
#else
	for(size_t r = 0; r < m_ObjNumber; r++)
#endif
	{
		if(!m_IsCenter[r])continue;
		for(size_t c = r; c < m_ObjNumber; c++)
		{
			if(!m_IsCenter[c])continue;
			if(r == c)
			{
				m_SimilarityMatrix[r][c] = DBL_MAX;
				continue;
			}
			m_SimilarityMatrix[r][c] = m_SimilarityMatrix[c][r] = EuDistanceSpecial(m_Centers[r], m_Centers[c]);
		}
	}
	delete []longLine;

#ifdef USE_EU_NO_COMBINE
Eu_Filter(part, distance);
return;
#endif

	m_TerminateDistance = distance;
	if(m_PartsNow < part)
	{
		cout << "Warning : " << " You're trying to partitioning " << m_PartsNow << " Objects into " 
				 << part << " parts!" << endl;
		cout << "Skip now..." << endl;
		part = m_PartsNow;
		return;
	}
	size_t toCombine1, toCombine2;
	while(1/*m_PartsNow > part*//* && minDis < m_TerminateDistance*/)
	{
		double minDis = DBL_MAX;
		// Find the most similar two ones.
#ifdef Force_Init_Seeds_For_Topo_Clutering
		for(size_t r = 1; r < m_ObjNumber; r++)
#else
		for(size_t r = 0; r < m_ObjNumber; r++)
#endif
		{
			if(!m_IsCenter[r])continue;
			for(size_t c = r; c < m_ObjNumber; c++)
			{
				if(!m_IsCenter[c])continue;
				if(m_SimilarityMatrix[r][c] < minDis)
				{
					minDis = m_SimilarityMatrix[r][c];
					toCombine1 = r;
					toCombine2 = c;
					break;
				}
			}
		}
		if(minDis > m_TerminateDistance) break;
		m_IsCenter[toCombine2] = false;

// Combine ..
		vector<Point> tmp;
		Point pnt1, pnt2;
		vector<double> len1, len2;
		computeLength(m_Centers[toCombine1], len1);
		computeLength(m_Centers[toCombine2], len2);
		for(size_t i = 0; i <= m_seg; i++)
		{
			getPointOnStreamline(m_Centers[toCombine1], len1, m_seg, i, pnt1);
			getPointOnStreamline(m_Centers[toCombine2], len2, m_seg, i, pnt2);
			tmp.push_back((pnt1 + pnt2)/2); 
		}
		m_Centers[toCombine1] = tmp;
// Update : m_SimilarityMatrix, partsNow, bigestDistance.
		for(size_t i = 0; i < m_ObjNumber; i++)
		{
			if(!m_IsCenter[i] || i == toCombine1)continue;
			double sim = 0;
			sim = EuDistanceSpecial(m_Centers[toCombine1], m_Centers[i]);
			m_SimilarityMatrix[i][toCombine1] = sim;
			m_SimilarityMatrix[toCombine1][i] = sim;
		}

		m_PartsNow--;
	}//while
  part = m_PartsNow;
}

vector<Point> *visHierarchicalClustering::getCenters()
{
	//size_t nonEmptyItem = 0;
	//for(size_t i = 0; i < m_ObjNumber; i++)
	//{
	//	if(m_IsCenter[i])nonEmptyItem++;
	//}
//	if(0 == nonEmptyItem);
	vector<Point> *tmpResult = new vector<Point>[m_PartsNow];
	size_t c = 0;
	for(size_t i = 0; i < m_ObjNumber; i++)
	{
		if(m_IsCenter[i])
		{
			// Find the most similar in m_streamlines...
			tmpResult[c++] = m_Centers[i];

//			double sim = DBL_MAX;
//			size_t assign = 0;
//			for(size_t j = 0; j < getStreamlinesNumber(); j++)
//			{
//				double tmp = EuDistance(j, m_Centers[i]);
//				if(tmp < sim)
//				{
//					sim = tmp;
//					assign = j;
//				}
//			}
//			for(size_t sz = 0; sz < getStreamlineSize(assign); sz++)
//			{
//				tmpResult[c].push_back(getPointOnStreamline(assign, sz));
//			}
//			c++;
		}
	}
	return tmpResult;
}

double visHierarchicalClustering::minDisSim(size_t objPre, size_t objNext)
{
	if( objPre == objNext ) return DBL_MAX;
	double *error = new double[m_Centers[objNext].size()];
	const double disError = pow(m_xmin-m_xmax,2)+pow(m_ymin-m_ymax,2)+pow(m_zmin-m_zmax,2);
	size_t count = 0;
	for(size_t i = 0; i < m_Centers[objNext].size(); i++)
	{
		nNearestNodes nNN(1);
		nNN.setDistError(disError);
		nNN.setSearchPnt( m_Centers[objNext][i].getX(),m_Centers[objNext][i].getY(),m_Centers[objNext][i].getZ() );
		lineTree[objPre].locateNodes(&nNN, 1);
		Point pntOnTree(nNN.index[1]->pos[0],nNN.index[1]->pos[1],nNN.index[1]->pos[2]);

		// direction
		size_t nearPointId = nNN.index[1]->id - 1;
		size_t preHead = nearPointId!=m_Centers[objPre].size()-1?nearPointId : m_Centers[objPre].size()-2;
		size_t preEnd = preHead + 1;
		size_t nextHead = i!=m_Centers[objNext].size()-1?i:m_Centers[objNext].size()-2;
		size_t nextEnd = nextHead + 1;
		Point direction1 = m_Centers[objPre][preEnd] - m_Centers[objPre][preHead];
		Point direction2 = m_Centers[objNext][nextEnd] - m_Centers[objNext][nextHead];

		double mod1 = pow( direction1.getX(), 2 ) + 
									pow( direction1.getY(), 2 ) +
									pow( direction1.getZ(), 2 );
		double mod2 = pow( direction2.getX(), 2 ) + 
									pow( direction2.getY(), 2 ) +
									pow( direction2.getZ(), 2 );
		mod1 = sqrt(mod1);
		mod2 = sqrt(mod2);
		if( mod1 < 1e-300 || mod2 < 1e-300 ) continue;
		direction1 /= mod1;
		direction2 /= mod2;

		double difference = direction1.getX() * direction2.getX() + 
												direction1.getY() * direction2.getY() +
												direction1.getZ() * direction2.getZ();
		if(difference < 0.1)
		{
			count++;
			if(10 <= count) break;
		}
		//
		error[i] = m_Centers[objNext][i].Distance(pntOnTree);
	}
	if(10 <= count)return DBL_MAX;
	double avg = 0, dist = 0;
	for(size_t i = 0; i < m_Centers[objNext].size(); i++)	{	avg += error[i]; }
	avg /= m_Centers[objNext].size();
//	for(size_t i = 0; i < m_Centers[objNext].size(); i++)	{	dist += abs(avg - error[i]); }
//	dist /= m_Centers[objNext].size();
	return sqrt(avg);
}

void visHierarchicalClustering::Eu_Filter(size_t &part, double distance)
{
	if(m_PartsNow < part) 
	{
		part = m_PartsNow;
		cout << "Ok!" << endl;
		return;
	}
	if(distance < 0) return;
	
	cout << "Filtering using eudistance no improvement..." << endl;
	double *accumDis = new double[m_ObjNumber];
	for(size_t i = 0; i < m_ObjNumber; i++) accumDis[i] = 0;

	for(size_t objPre = 0; objPre < m_ObjNumber; objPre++)
	{
		if( !m_IsCenter[objPre] ) continue;
		for(size_t objNext = objPre + 1; objNext < m_ObjNumber; objNext++)
		{
			if(   !m_IsCenter[objNext]
			    || objPre == objNext 
//					|| accumDis[objNext] > 2 * distance
			/*
			 *  || m_Length[objPre] < m_Length[objNext]
			 *	|| objNext == 0 
			*/
				) continue;
			double dist = EuDistance( m_Centers[objPre], m_Centers[objNext] );	
			if(dist < distance)
			{
				m_IsCenter[objNext] = false;
				m_PartsNow--;
				accumDis[objPre] += dist;
				if( m_PartsNow <= part )
				{
					part = m_PartsNow;
					cout << "ok ,return..." << endl;
					return;
				}
			}
		}
	}

	part = m_PartsNow;
	delete []accumDis;
	cout << "Normally return..." << endl;
}

double visHierarchicalClustering::EuDistanceSpecial(strmLine &line1, strmLine &line2)
{
	double sim = 0;

	double *error = new double[m_seg + 1];
	size_t count = 0;
	for(unsigned int i = 0; i <= m_seg; i++)
	{
		error[i] = line1[i].Distance(line2[i]);
//Compute the direction difference...
/*    	line1  preHead ---> preEnd, line2  nextHead ---> nextEnd */
		size_t head = (m_seg == i)?(m_seg-1):i;
		size_t end = head + 1;
		Point direction1 = line1[end] - line1[head];
		double mod1 = pow( direction1.getX(), 2 ) + 
									pow( direction1.getY(), 2 ) +
									pow( direction1.getZ(), 2 );
		Point direction2 = line2[end] - line2[head];
		double mod2 = pow( direction2.getX(), 2 ) + 
									pow( direction2.getY(), 2 ) +
									pow( direction2.getZ(), 2 );
		mod1 = sqrt(mod1);
		mod2 = sqrt(mod2);
		if( mod1 < 1e-300 || mod2 < 1e-300 ) continue;
		direction1 /= mod1;
		direction2 /= mod2;
		double difference = direction1.getX() * direction2.getX() + 
												direction1.getY() * direction2.getY() +
												direction1.getZ() * direction2.getZ();
		if(difference < 0.5)
		{
			count++;
			if(5 <= count) return DBL_MAX;
		}
	}
	double avg = 0, dist = 0;
	for(size_t i = 0; i <= m_seg; i++) {	avg += error[i]; }
	avg /= (1 + m_seg);
	for(size_t i = 0; i <= m_seg; i++) { dist += abs( avg - error[i] ); }
	dist /= (1 + m_seg);

	sim = 0.65 * dist + 0.35 * avg;
	delete []error;
	return sim;
}

void visHierarchicalClustering::filter(size_t &part, double distance)
{
	if(m_PartsNow < part) 
	{
		part = m_PartsNow;
		cout << "Ok!" << endl;
		return;
	}
	if(distance < 0) return;
	
	buildMinDisTree();
	cout << "Filtering..." << endl;
	const double disError = pow(m_xmin-m_xmax,2)+pow(m_ymin-m_ymax,2)+pow(m_zmin-m_zmax,2);
	double *accumDis = new double[m_ObjNumber];
	for(size_t i = 0; i < m_ObjNumber; i++) accumDis[i] = 0;

#ifdef Force_Init_Seeds_For_Topo_Clutering
	for(size_t objPre = 1; objPre < m_ObjNumber; objPre++)
#else
	for(size_t objPre = 0; objPre < m_ObjNumber; objPre++)
#endif
	{
		if( !m_IsCenter[objPre] ) continue;
#ifdef Force_Init_Seeds_For_Topo_Clutering
		for(size_t objNext = 1; objNext < m_ObjNumber; objNext++)
#else
		for(size_t objNext = 0; objNext < m_ObjNumber; objNext++)
#endif
		{
			if(   !m_IsCenter[objNext]
			    || objPre == objNext 
//					|| accumDis[objNext] > 5 * distance
			/*
			 *  || m_Length[objPre] < m_Length[objNext]
			 *	|| objNext == 0 
			*/
				) continue;
			double *error = new double[m_Centers[objNext].size()];
			size_t count = 0;
			for(size_t i = 0; i < m_Centers[objNext].size(); i++)
			{
				nNearestNodes nNN(1);
				nNN.setDistError(disError);
				nNN.setSearchPnt( m_Centers[objNext][i].getX(),m_Centers[objNext][i].getY(),m_Centers[objNext][i].getZ() );
				lineTree[objPre].locateNodes(&nNN, 1);
				Point pntOnTree(nNN.index[1]->pos[0],nNN.index[1]->pos[1],nNN.index[1]->pos[2]);
//Compute the direction difference...
/*    	line1  preHead ---> preEnd, line2  nextHead ---> nextEnd */
				size_t nearPointId = nNN.index[1]->id - 1;
				size_t preHead = nearPointId!=m_Centers[objPre].size()-1?nearPointId : m_Centers[objPre].size()-2;
				size_t preEnd = preHead + 1;
				size_t nextHead = i!=m_Centers[objNext].size()-1?i:m_Centers[objNext].size()-2;
				size_t nextEnd = nextHead + 1;
				Point direction1 = m_Centers[objPre][preEnd] - m_Centers[objPre][preHead];
				double mod1 = pow( direction1.getX(), 2 ) + 
											pow( direction1.getY(), 2 ) +
											pow( direction1.getZ(), 2 );
				Point direction2 = m_Centers[objNext][nextEnd] - m_Centers[objNext][nextHead];
				double mod2 = pow( direction2.getX(), 2 ) + 
											pow( direction2.getY(), 2 ) +
											pow( direction2.getZ(), 2 );
				mod1 = sqrt(mod1);
				mod2 = sqrt(mod2);
				if( mod1 < 1e-300 || mod2 < 1e-300 ) continue;
				direction1 /= mod1;
				direction2 /= mod2;
				double difference = direction1.getX() * direction2.getX() + 
														direction1.getY() * direction2.getY() +
														direction1.getZ() * direction2.getZ();
				if(difference < 0.1)
				{
					count++;
					if(10 <= count) break;
				}
//Compute the direction difference

				error[i] = m_Centers[objNext][i].Distance(pntOnTree);
			}
			if(count >= 10)continue;
			double avg = 0, dist = 0;
			for(size_t i = 0; i < m_Centers[objNext].size(); i++)	{	avg += error[i]; }
			avg /= m_Centers[objNext].size();
			for(size_t i = 0; i < m_Centers[objNext].size(); i++)	{	dist += abs(avg - error[i]); }
			dist /= m_Centers[objNext].size();
			dist =  0.3*dist + 0.7*avg;

			if(dist < distance)
			{
				m_IsCenter[objNext] = false;
				m_PartsNow--;
				accumDis[objPre] += dist;
				if( m_PartsNow <= part )
				{
					part = m_PartsNow;
					cout << "ok ,return..." << endl;
					return;
				}
			}
			delete []error;
		}
	}

	part = m_PartsNow;
	delete []accumDis;
	cout << "Normally return..." << endl;
}

void visHierarchicalClustering::buildMinDisTree()
{
	lineTree = new kdTree[m_ObjNumber];
	cout << "Building the tree..." << endl;
	for(size_t i = 0; i < m_ObjNumber; i++)
	{
		if(!m_IsCenter[i])continue;
		lineTree[i].set(m_Centers[i].size());
		for(size_t sz = 0; sz < m_Centers[i].size(); sz++)
		{
			lineTree[i].store( m_Centers[i][sz].getX(), 
												 m_Centers[i][sz].getY(), 
												 m_Centers[i][sz].getZ() 
											 );
		}
		lineTree[i].treeBalance();
	}
}

void visHierarchicalClustering::cylinder()
{
	for(size_t i = 1; i <= m_dataField->getPntNumberOfModel(); i++)
	{
		vector<double> coords, dataV;
		m_dataField->getPointCoordinates(i, coords);			
		m_dataField->getDataValueOfPoint(i, dataV);
		if( abs(dataV[0]) < 1e-6 && abs(dataV[0]) < 1e-6 && abs(dataV[0]) < 1e-6 )
		{
			Point tmp(coords[0], coords[1], coords[2]);
			cords.push_back(tmp);
		}
	}
}

}

